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We solve the time-dependent Schrodinger equation for the combination of a spin system 
interacting with a spin bath environment. In particular, we focus on the time development 
of the reduced density matrix of the spin system. Under normal circumstances we show 
that the environment drives the reduced density matrix to a fully decoherent state, and 
furthermore the diagonal elements of the reduced density matrix approach those expected 
for the system in the canonical ensemble. We show one exception to the normal case is if 
the spin system cannot exchange energy with the spin bath. Our demonstration docs not 
rely on time-averaging of observablcs nor does it assume that the coupling between system 
and bath is weak. Our findings show that the canonical ensemble is a state that may result 
from pure quantum dynamics, suggesting that quantum mechanics may be regarded as the 
foundation of quantum statistical mechanics. 
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1. Introduction 

Statistical mechanics is one of cornerstones of modern physics but its foundations and basic 
postulates are still under debate. ^"^^^ There is a common believe that a generic "system" that 
interacts with a generic environment evolves into a state described by the canonical ensemble. 
Experience shows that this is true but a detailed understanding of this process, which is 
crucial for a rigorous justification of statistical physics and thermodynamics, is still lacking. 
In particular, in this context the meaning of "generic" is not clear. The key question is to 
what extent the evolution to the equilibrium state depends on the details of the dynamics of 
the whole system. 

Earlier demonstrations that the system can be in the canonical ensemble state are based 
on showing that time-averages of the expectation dynamical variables of the system approach 
their values for the subsystem that is the thermal equilibrium state^"^^ or do not consider the 
dynamics of the system but assume that the state of the whole system has a special property 
called "canonical typicality" ^"^■^^ in which case it is as yet unclear under which conditions the 
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whole system will evolve to the region in Hilbert space where its subsystems are in the thermal 
equilibrium state. A very different setting to study nonequilibrium quantum dynamics is to 
start from an eigenstate of some initial Hamiltonian and push the system out of this state by a 
sudden change of the model parameters.^^"^^^ To the best of our knowledge, it has not yet been 
shown that this approach leads to the establishment of the canonical equilibrium distribution. 
Finally, we want to draw attention to the fact that a demonstration of relaxation to the 
canonical distribution requires a system with at least three different eigenenergies because 
a diagonal density matrix of a two-level system can always be represented as a canonical 
distribution.^'^' 

The main result of this paper is that we show, without any time-averaging procedure 
or any approximation, that systems embedded in a closed quantum system generally evolve 
to their canonical distribution states. This result complies with the fact that if we make a 
real measurement of a thermodynamic property, we observe its equilibrium value without 
having to perform time averaging. Furthermore, we show that the relaxation to the canonical 
distribution is not limited to the regime of weak coupling between system and environment, 
an assumption that is often used.^'^~^^^ 

2. General theory 



In general, the state of a closed quantum system is described by a density matrix. ^^'^^^ 
The canonical ensemble is characterized by a density matrix that is diagonal with respect to 
the eigenstates of the system Hamiltonian, the diagonal elements taking the form exp{—f3Ei) 
where (3 = l/ksT is proportional to the inverse temperature (ks is Boltzmann's constant) and 
the Ei^s denote the eigenenergies. The time evolution of a closed quantum system is governed 
by the time-dependent Schrodinger equation (TDSE).^^'^^) If the initial density matrix of an 
isolated quantum system is non-diagonal, then, according to the TDSE, its density matrix 
remains nondiagonal and never approaches the thermal equilibrium state with the canonical 
distribution. Therefore, in order to thermalize the system S, it is necessary to have the system 
S interact with an environment {E), also called heat bath. Thus, the Hamiltonian of the whole 
system {S + E) takes the form H = Hg + He + Hse, where Hs and He are the system and 
environment Hamiltonian, respectively and Hse describes the interaction between the system 
and environment. 

The state of system S is described by the reduced density matrix 

p{t) = T^Ep{t), (1) 

where p {€) is the density matrix of the whole system at time t and Tr^ denotes the trace over 
the degrees of freedom of the environment. The system S is in its thermal equilibrium state 
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if the reduced density matrix takes the form 

p=e-^^^/Tr5e-/^^^ (2) 

where Trs denotes the trace over the degrees of freedom of the system S. Therefore, in order 
to demonstrate that the system S, evolving in time according to the TDSE, relaxes to its 
thermal equilibrium state one has to show that p{t) p for t > to where to is some finite 
time. 

The difference between the state p{t) and the canonical distribution jo is most conveniently 
characterized by the two quantities S{t) and a{t) defined by 
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Here N denotes the dimension of the Hilbert space of system S and Pij{t) is the matrix 
element of the reduced density matrix pin the representation that diagonalizes Hg- As 
the system relaxes to its canonical distribution both S{t) and a{t) vanish, b{t) converging to 
p. As a{t) is a global measure for the size of the off-diagonal terms of the reduced density 
matrix, a{t) also characterizes the degree of coherence in the system: If a{t) =0 the system 
is in a state of full decoherence. 

3. Model and simulation method 

To study the evolution to the canonical ensemble state in detail, we consider a general 
quantum spin- 1/2 model defined by the Hamiltonians 

n s — 1 lis 

Hs = E E J^srsf, (6) 

i=l j=i+la=x.y,z 
n—1 n 

He = -EE H ^Jirir^ (7) 

i=l j=i+l a=x,y,z 

ns n 

= -EE E ^.spf. (8) 

i=l j=l a=x,y,z 

Here the S^'s and /"'s denote the spin- 1/2 operators of the system and environment re- 
spectively (we use units such that h and ks are one). Analytic expressions for p{t) can only 
be obtained for very special choices of the exchange integrals J"^, and A"^. but it is 
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straightforward to solve the TDSE numerically for any choice of the model parameters. Here, 
we numerically solve the TDSE for H = Hg + He + Hse using the Chebyshev polynomial 
algorithm. These ab initio simulations yield results that are very accurate (at least 10 
digits), independent of the time step used.^^^ 

The state, that is the density matrix p{t) of the whole system at time t is completely 
determined by the choice of the initial state of the whole system and the numerical solution 
of the TDSE. In our work, the initial state of the whole system (S+E) is a pure state. This 
state evolves in time according to 

Mt)) = e-^*|*(0)) = ^f^c(i,p,t)K,p), (9) 

i=i p=i 

where the states {\i,p)} denote a complete set of orthonormal states. In terms of the expansion 
coefficients c{i,p,t), the reduced density matrix reads 

2" 2" 

p(t)ij = TrE'^'^c*{i,q,t)c(j,p,t)\j,p){i,q\ 

p=l q=l 

2" 

= J^c*(i,p,t)c(j,p,t), (10) 

p=i 

which is easy to compute from the solution of the TDSE. Another quantity of interest that 
can be extracted from the solution of the TDSE is the local density of states (LDOS) 

1 r+oo 

LDOS(£:) = — / dt e-'^\^{0)\e-'"^\^{0)) 

fc=i 

2"s 2" 

= ^i;c*(z,|.,0)c(i,^>,t), (11) 
1=1 p=l 

where D = 2"'+"s'^ {Iv'fc)}) and {-Efe} denote the dimension of the Hilbert space, the eigenstates 
and eigenvalues of the whole system, respectively. The LDOS is "local" with respect to the 
initial state: It provides information about the overlap of the initial state and the eigenstates 
oiH. 

The notation to specify the initial state is as follows: \GROU N D) g is the ground state or 
a random superposition of all degenerated ground states of the system; \RANDOM) g denotes 
a random superposition of all possible basis states; \UU)g is a state in which all spins of the 
system are up meaning that in this state, the expectations value of each spin is one; \UD)g 
is a state in which two nearest-neighbor spins of the system are antiparallel implying that in 
this state, the correlation of their z-components is minus one; and \RR)s denotes the product 
state of random superpositions of the states of the individual spins of the system. The same 
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Fig. 1. (Color online) Simulation results for the diagonal elements pi = pii{t) of the density matrix 
of S, the energy Es = Es{t), the effective inverse temperature b = b{t) and its variance S = S{t), 
and a = a{t) which is measure for the decoherence in S, as obtained by solving the TDSE for the 
whole system with Heisenberg-ring Hs (J = —1, ns = 4), Heisenberg-type Hse (A = 0.3), spin 
glass He ^ 1, n = 18), and r = tt/IO. The initial state of the whole system is a product state 
of \UD)g and \RANDOM)j^. 



notation is used for the spins in the environment, the subscript S being replaced by E. 

As we report results for many different types of spin systems it is useful to introduce 
a simple terminology to classify them according to symmetry and connectivity. The terms 
"XY" , "Heisenberg" , "Heisenberg-type" and "Ising" system refer to the cases J^j = ■ = J 
and Jfj- = 0, Jfj = Jfj = J^j = J, Jij uniform random in the range [— |^|,|^|], and 
J[j = Jfj = and Jfj = J, respectively. The same terminology of symmetry is used for 
the Hamiltonian He of the environment and for the interaction Hamiltonian Hse- In our 
model, all the spins of the system interact with each spin of the environment. To characterize 
the connectivity of spins within the system (environment), we use the term "ring" for spins 
forming a one-dimensional chain with nearest-neighbor interactions and periodic boundary 
conditions, "triangular-lattice" if the spins are located on a two-dimensional triangular lattice 
with nearest-neighbor interactions, and "maximum-connectivity-system" when all the spins 
within the system (environment) interact with each other. 
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Fig. 2. (Color online) Simulation results for the energy Es = Es{t), the effective inverse temperature 
b = b(t), its variance S = S{t), and the deviation from a diagonal matrix cr = a{t) as obtained 
by the solution of the TDSE for a variety of different systems S coupled to a spin glass He via 
a Heisenberg-type Hse- The systems used are a: XY-ring, b and f: Heisenberg-ring, c: Ising- 
ring, d: Heisenberg-triangular-lattice, and e: spin glass (Heisenberg-type maximum-connectivity- 
system). The initial states of the whole system are a: \GROUND)s ® \RANDOM) e, b: \UD)s ® 
\RANDOM)e, c: \UU)s ® \RR)e, d: \UU)s \RANDOM)e, e: \GROUND)s ® \UD)e, and 
f: \UD)s <8 \GROU N D) E- The numbers of spins in the system are rig = 8 for cases a-c and 
ns — Q for cases d-f. The numbers of spins in the environment is n = 16 for all cases. The model 
parameters are J = —1, A = 0.3 and = 1, except for case e in which A = 1. 

4. Results 

In earlier work, it was found that a frustrated spin glass (Heisenberg-type maximum- 
connectivity-system) environment is very effective for creating full decoherence {a — > 0) in 
a two-spin system. As o" — > is a necessary condition for the state of the system to 
converge to its canonical distribution, we have chosen spin glass environments, which have no 
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obvious symmetries, for further exploration. 

First, we consider a system {Hs- Heisenberg-ring) interacting (Hse- Heisenberg-type) 
with an environment (He- spin glass). The system has four distinct eigenvalues (Ei = —2, 
£^2-4 = —1) -£'5-11 = 0, and £'i2_i6 = 1) and sixteen different eigenstates. The environment 
has 2^^ eigenstates. During the time-integration of the TDSE, the reduced density matrix of 
the system is calculated every r = vr/lO. Following the general procedure described earlier, the 
values of the diagonal elements pa yield an estimate for the effective inverse temperature b{t), 
the error S{t) for this estimate and the measure a{t) for the deviation from a non-diagonal 
matrix. We also monitor the energy Es{t) = Trsp{t)Hs, of the system. 

From the simulation results, shown in Fig. 1, it is clear that for t > 50r, each diagonal 
element pa of the reduced density matrix converges to one out of four stationary values, 
corresponding to the four non-degenerate energy levels of the system. This convergence is a 
two-step process. First the system looses all coherence, as indicated by the vanishing of a (t) 
for t > 50r. The time dependence of a {t) fits very well to an exponential law 

<T (t) = (Too + Ae-*/^^ (12) 

with (Too = 0.00128, A = 0.602 and T2 = S.Olr. Likewise, the vanishing of S{t) on the 
same time-scale indicates that the density matrix of the system converges to the canonical 
distribution. The effective temperature b{t) and the energy of the system Es {t) also fit very 
well to the exponential laws 

h{t) = l3 + Be-''^\ (13) 

and 

E{t) = E^ + Ce-'/'^\ (14) 

with 13 = 0.0962, B = -0.900, and Ti = 13. 3t and E^^ = -0.0745, C = -0.952. The 
estimated values for Ti and T2 change very little if we choose different random realizations 
for the initial state of the environment or for the model parameters and A"^- (data not 
shown) but if we change their range, Ti and T2 also change, as naively expected. 

The simulation demonstrates that the system first looses all coherence and then, on a 
longer time-scale, relaxes to its thermal equilibrium state with a finite temperature. In terms 
of the theory of magnetic resonance,^^) Ti and T2 axe the times of dissipation and dephasing, 
respectively. Note that in contrast to the cases considered in the theory of nuclear magnetic 
resonance, in most of our simulations, Hs, He and Hse are comparable so the standard 
perturbation derivation of a and E does not work. In the case of very small He, one should 
expect, instead of an exponential decay of a and E, a Gaussian decay, as observed in our 
earlier work.^^~^^) 

Results for systems [Hs) with different symmetries and connectivities that interaction 
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Fig. 3. (Color online) Simulation results for the local density of states as a function of the energy. 
Solid line: Case corresponding to Fig. 2b. The initial state is \UD)s ® \RANDOM) e; Dashed 
line: Case corresponding to Fig. 2f. The initial state is \UD)s <8) \GROUND)e- 



with the same type of environments {He) via the same type of couphngs {Use) are shown 
in Fig. 2. The systems used are an XY-ring, a Heisenberg-ring, an Ising-ring, a Heisenberg- 
triangular-lattice, and a spin glass. From Fig. 2, it is clear that independent the internal 
symmetries and connectivity of the system and independent the initial state of the whole 
system (except for case f in which the environment is initially in its ground state), all systems 
relax to a state with full decoherence. Notice that in case b, a vanishes exponentially with 
time, whereas in other cases (a,c,d,e), a initially increases and then vanishes exponentially 
with time, due to the entanglement between the system and the environment. This observation 
is in concert with our earlier work.^^^^^^ 

Furthermore, in all cases except f, the system always relaxes to a canonical distribution 
{5 — > 0) as soon as it is in the state with full decoherence {a — > 0), indicating that the 
time of decoherence {T2) and the thermalization is almost the same. In agreement with the 
results depicted in Fig. 1, the decoherence time T2 is shorter than the typical time scale Ti 
on which the system and environment exchange energy and the effective inverse temperature 
h{t) reaches its stationary value. 

The case f is easily understood in terms of the local density of states. In Fig. 3 we show 
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the LDOS for the cases b and f, the only difference between these two cases being the initial 
state of the environment. Up to a trivial normalization factor, the LDOS curve for case b 
is indistinghuisable from the density of states (data not shown) calculated from the solution 
of the TDSE using the technique described in Ref.^^^ This suggests that if the environment 
starts from the random superposition of all its states, all states of the whole system may 
participate in the decoherence/relaxation process. In contrast, the LDOS curve for case f has 
a very small overlap with the density of states (the curve of which coincides with the solid 
line in Fig. 3). Therefore, starting with an environment in the ground state, only a relatively 
small number states participates in the decoherence process, as confirmed by the results for 
a{t) shown in Fig. 2f. 

For completeness, we discuss a two other situations in which, for fairly obvious reasons, 
the system cannot relax to its canonical distribution. Obviously, if the energy of the system 
is conserved {[Hs, H] =0), the system cannot exchange energy with the environment and we 
should not expect relaxation to the canonical distribution. In this case, as shown in Fig. 4, 
after the system S has reached a state with full decoherence, its density matrix does not 
converge to the canonical state. Likewise, if the range of energies of the environment E is 
too small compared to that of the system <^ | J[) as in the example shown in Fig. 4b, 
there is no convergence to the canonical state either. It is to be noted that in both cases, the 
interaction with the environment leads to perfect decoherence {cr{t) ~ 0, see insets) such that 
the reduced density matrix converges to a diagonal matrix. However, from Fig. 4, it follows 
that S relaxes to a kind of microcanical state in which the states in each energy subspace have 
equal probability, the probabilities to end up in a subspace depending on the initial state. 

Disregarding the three cases mentioned earlier, the simulation results presented in Figs. (1) 
and (2) suggest that the state of a system generally relaxes to the canonical distribution when 
the system is coupled to an environment of which the dynamics is sufficiently complex also 
in the case that the interaction between system and environment cannot be regarded as 
a perturbation. There are exceptions but these are easily understood: Either there are not 
enough states available for the decoherence (Fig. 2f ) to yield a diagonal reduced density matrix 
or the energy relaxation (Fig. 4) is not effective in letting the diagonal reduced density matrix 
relax to the canonical distribution. 

Although we have only presented results for a spin glass environment He, our results (not 
shown) for any of the choices for Hs and He mentioned earlier, in combination a Heisenberg- 
type HsE interaction between system and environment, or for Hs and Hse in combination 
Heisenberg-type He leads to the same conclusion, namely that the state of a system relaxes 
the canonical distribution. 
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5. Discussion 

The results presented here have been obtained from an ab initio numerical solution of 
the TDSE in the absenee of, for instance, dissipative mechanisms, and demonstrate that the 
existence of the canonical distribution, a basic postulate of statistical mechanics, is a direct 
consequence of quantum dynamics. 

We have shown that if we have a system S that interacts with an environment E and the 
whole system S + E forms a closed quantum system that evolves in time according to the 
TDSE, S and E can exchange energy, the range of energies of E is large compared to the 
range of energies of S, and the interaction between S and E leads to full decoherence of S, 
then the state of S relaxes to the canonical distribution. Note that only the condition of full 
decoherence is a nontrivial requirement. 

We emphasize that our conclusion does not rely on time averaging of observables, in concert 
with the fact that real measurements of thermodynamic properties yield instantaneous, not 
time-averaged, values. Furthermore and perhaps a little counter intuitive, our results show 
that relatively small environments (~ 20 spins) are sufficient to drive the system S to thermal 
equilibrium and that there is no need to assume that the interaction between the system and 
environment is weak, as is usually done in kinetic theory. 

In conclusion: The work presented here strongly suggests that the canonical ensemble, 
being one of the basic postulates of statistical mechanics, is a natural consequence of the 
dynamical evolution of a quantum system. This conclusion may be exciting but as quantum 
mechanics describes the dynamics of a system and statistical mechanics gives us the distribu- 
tion when the system is in the equilibrium state, these two successful theories should not be 
in conflict once the conditions for the system to relax to its thermal equilibrium are satisfied. 

6. Aknowledgments 

It is a pleasure to thank S. Miyashita, F. Jin, S. Zhao and K. Michielsen for many helpful 
discussions. We are grateful to S. Miyashita and M. Novotny for several suggestions to improve 
the manuscript. This work was partially supported by NCF, The Netherlands. 



10/?? 



J. Phys. Soc. Jpn. 



Full Paper 




20 40 60 80 100 




50 100 150 200 250 

Vt 

Fig. 4. (Color online) Simulation results for a Heisenberg-ring Hg (J = —5, ng = 4, initial state 
\UD)g) coupled to a spin glass He — 0.15, n — 16, initial state \RANDOM) e) via (a) 
Heisenberg HgE (A = 0.075) or (b) Heisenberg-type Hse (A — 0.15). Although full decoherence 
is observed in both cases, the the system S only relaxes to a state with equal probabilities within 
each energy subspace, that is to a "microcanonical" state per energy subspace. 
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